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Abstract 

We present results of Brownian dynamics simulations of tethered nanospheres and tethered 
nanorods. Immiscibility between tether and nanoparticle facilitates microphase separation into 
the bicontinuous, double gyroid structure (first reported by Iacovella et al. [Phys. Rev. E 75 
(2007)] and Horsch et al [J. Chem. Phys. 125 (2006)] respectively). We demonstrate the ability 
of these nanoparticles to adopt distinct, minimal energy local packings, in which nanospheres form 
icosahedral-like clusters and nanorods form splayed hexagonal bundles. These local structures 
reduce packing frustration within the nodes of the double gyroid. We argue that the ability to 
locally order into stable structures is key to the formation of the double gyroid phase in these 
systems. 



PACS numbers: 



Block copolymers and surfactants have long been known to self- assemble into a wide vari- 
ety of complex structures where the assembly is driven by immiscibility between chemically 
distinct blocks in the polymers [I] and between distinct head and tail groups in surfactants 
[2]. These ordered structures are highly sought for applications at the nanoscale, ranging 
from photonic-bandgap materials [3j to templates for nanoparticle assembly [4] and hydrogen 
storage. Hybrid building blocks have recently been created that resemble block copolymers 
where the individual blocks consist of nanoparticles and polymers j5j EH d El |9]. These 
hybrid building blocks, or tethered nanoparticles, constitute a class of "shape amphiphiles" 
[TUl [TT] where microphase separation occurs due to the immiscibility between the nanopar- 
ticle and polymeric tether resulting in mesostructured equilibrium phases that resemble the 
morphologies of block copolymers and surfactants [T21 [13] . 

In previous work, we examined the interplay between microphase separation and nanopar- 
ticle geometry for tethered nanospheres (TNS) [TSKH] and tethered nanorods (TNR) [T2l[T5] . 
finding unique changes to phase behavior as compared to flexible surfactants and block 
copolymers. For example, the phase behavior of the nanosphere- aggregating TNS system 
includes hexagonally packed cylinders, the double gyroid morphology, and perforated lamel- 
lar phases where there is a predominance towards icosahedral ordering of nanospheres [14] 
and a lamellar phase with HCP ordering of nanospheres [HJ. The tether- aggregating TNS 
system does not form the double gyroid structure, as would be expected for surfactants, 
instead forming perforated lamella [13]. The phase behavior of the nanorod-aggregating 
TNR system includes a hexagonally packed cylinder phase where the nanorods twist along 
the length of the cylinder, hexagonally and tetragonally perforated lamellar phases where 
the nanorods form a smectic structure, and a lamellar phase with smectic C-like ordering 
of the nanorods [12]. A transition from hexagonally packed cylinders to the double gyroid 
morphology is seen upon reducing the length of the tether in the TNR system [15J, similar 
to what was observered for rod-coil liquid crystals [16J. 

In this work, we examine and compare the double gyroid (DG) microstructure formed 
by the TNS and TNR building blocks with attractive nanoparticles and repulsive tethers in 
order to learn about the stability of the DG structure. The DG is a bicontinous structure 
where the nanoparticles form two distinct, interpenetrating networks. The DG structure is 
of particular interest as it is seen as a candidate for catalytic materials, high conductivity 
nanocomposites [IT], and photonics applications [3]. In this work we first present our TNS 
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and TNR models and simulation method in section [T| followed by the calculation of the 
Flory-Huggins interaction parameter to allow for comparisons between the TNS and TNR 
systems in section [TTJ. In section |III A| we investigate packing frustration within the DG 
and its connection to particle geometry. In sections |III B| and HI C we investigate the local 



configurations of the nanospheres and nanorods respectively. In section |IV| we provide 
concluding remarks. 



I. MODEL AND METHOD 



To study tethered nanoparticles, we consider a general class of tethered nanoparticles 
rather than any one specific system and use empirical pair potentials that have been suc- 
cessful in the study of block copolymers and surfactants [18]. We utilize minimal models 
that capture the essential physics of the problem, specifically the geometry of the nanoparti- 
cle, immiscibility between tether and nanoparticle, and flexibility of the polymer tether. We 
examine tethered nanospheres and tethered nanorods in selective solvent and additionally 
explore the connection between these building blocks and diblock copolymers in selective 
solvent. The natural units of these systems are: a, the diameter of a tether bead; m, the 
mass of a tether bead; and e, the Lennard- Jones well depth. Bulk system volume fraction, 

is defined as the ratio of volume of the beads to the system volume, the dimensionless 
time is i*=a^/m/e, and the degree of immiscibility and solvent quality are determined by 
the inverse temperature, 1/T* — e/kBT. 



A. Tethered Nanospheres 

Nanospheres are modeled as beads of diameter 2.0a connected to tethers via finitely 
extensible non-linear elastic (FENE) springs [19J. Tethers are modeled as bead-spring chains 
containing eight beads of diameter a connected via FENE springs; a schematic of the model 
building block is shown in Figure [T] To model the attractive interaction between NPs we 
use the Lennard- Jones potential (LJ) where particle-particle interactions are shifted to the 
surface (Equation [I]); for nanoparticle-nanoparticle interactions We Set T shift — (2.0a - a) 
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FIG. 1: Model building blocks utilized. 

and r cutoff = 2.5a + r shift . 
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, r > r ^toff 

Solvophilic tethers interact via the purely repulsive Weeks-Chandler- Andersen (WCA) soft- 
sphere potential to account for short-range, excluded volume interactions (Equation [2]); for 
tether-tether interactions r s hift = and r cu toff = 2 1/6 cr. 



, r > r cuto// 

Nanosphere-tether interactions are treated with the purely repulsive WCA soft-sphere poten- 
tial to account for short-range, excluded volume interactions (Equation [2]); for nanosphere- 
tether interactions r shift — §(2.0cr — a) and r cuto ff = 2 1 / 6 a + r shift- This model relates well 
to experimentally synthesized building blocks including tethered quantum dots [5], polymer- 
functionalized fullerenes [6], tethered nanospheres formed by crosslinking one block of a BCP 
[7j, and divalent nanospheres [9]. 
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B. Tethered Nanorods 



Nanorods are modeled as a rigid collection of five beads of diameter a connected to tethers 
via FENE springs; a schematic of the model building block is shown in Figure [TJ The beads 
in the rod are spaced a distance a apart creating a relatively "rough" rod for computational 
efficiency. In previous studies the assembly behavior of a "rough" rod was indistinguishable 
from that of a "smooth" rod [15]. Moreover, this rough rod relates well to colloidal rod 
assemblies that have recently been fabricated [20j [21] . Tethers are modeled as bead-spring 
chains containing two beads of diameter a connected via FENE springs. Interactions between 
nanorods are modeled by the LJ potential (Equation [I]), where r s hift = and r cuto ff = 2.5a. 
Solvophilic tethers and species of different type interact via the WCA potential (Equation 
[2]) to account for short-range, excluded volume interactions, where in all cases r s hift = and 
r cutoff = 2 1/6 a. 

C. Block Copolymers 

Diblock copolymer systems (BCPs) are modeled as bead spring chains where individual 
beads of diameter a are connected via FENE springs; a schematic of the model building block 
is shown in Figure [l] Particles in the A-block (analogous to the head group in the TNS 
and TNR systems) are modeled as five beads that interact via the LJ potential (Equation 
[T]), where r shift = and r cuto ff = 2.5a. Particles in the B-block (analogous to the tethers 
in the TNS and TNR systems) are modeled by two beads of diameter a. Solvophilic B- 
blocks (tethers) and species of different type interact via the WCA potential (Equation [2]) 
to account for short-range, excluded volume interactions where in all cases r s hift = and 
rcutoff = 2 1/6 a. 

Additional explanation and utilization of these models can be found in references fTTl [T4] 
for TNS, [IU H21 H5] for TNR, and [EH [22] for BCP. 

D. Method of Brownian Dynamics 

To realize the long time scales and large systems required to study the self-assembly of 
complex mesophases, we use the method of Brownian dynamics (BD). The trajectory of 
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each bead is governed by the Langevin equation: 



mMt) = F?{n(t)) + Ff (t) - jVi(t) 



(3) 



where m i: r i: v$, Ff , Ff and 7^ correspond to the mass, position, velocity, conservative force, 
random force, and friction coefficient of bead z, respectively. We assume that there are no 
spatial or temporal fluctuations in the friction coefficient and fix 7^ = 1.0, which limits the 
range of ballistic motion of a bead to approximately 1.0a. The random force is independent 
of the conservative force and satisfies the fluctuation-dissipation theorem: 



The friction coefficient and random force act as a non-momentum-conserving heat bath; 
the combination of these two terms helps to minimize numerical roundoff errors that can 
occur over long simulations runs. The stationary solution of the Langevin equation is the 
Boltzmann distribution and therefore BD samples the canonical (NVT) ensemble. We do 
not explicitly include solvent particles; however, the frictional and random forces help to 
implicitly account for some of the effects of solvent. To incorporate rotational degrees 
of freedom into nanorods, we utilize the equations of rotation for linear bodies [23j. In 
all cases, particle beads are advanced through time using the leapfrog algorithm with a 
timestep At = 0.01. The general procedure employed is to start initially at an athermal 
condition, where all interactions are treated using the WCA potential, and allow the system 
to become well mixed. For a fixed 0, we then incrementally cool the system towards a 
final target temperature T* with selective solvent interactions, allowing it to run for several 
million timesteps at each intermediate T*. We also perform heating runs, where we start 
from the final structure and incrementally raise T* until we reach a disordered structure. 
To help avoid kinetically arrested structures, we repeat the simulations, varying the cooling 
sequence; for each building block, three unique cooling sequences were used. Approximately 
50 individual points were simulated for the TNS DG phase and 40 for the TNR DG phase, 
in total using approximately 10000 cpu hours. Simulations were performed with 500 TNS 
building blocks and 800 TNR building blocks, as these correspond to the unit cell of the DG 
morphology. 



<F?(f)> = 

{Ff(t)Ff(t'))=6 7 k B TS i ,S(t-t') 



(4) 
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II. FLORY-HUGGINS X PARAMETER 



Although the models for the TNS, TNR and BCP are very similar, key parameters such as 
T* do not necessarily correspond to the same statepoint. For example, in Reference [12] the 
rigidity of the tethered nanorod induced phase separation at a higher T* than the equivalent 
coil-coil BCP. Additionally, it is well known that the order-disorder temperature in BCP 
systems will scale as the number of beads [21] , so we expect that a tethered nanorod system 
with a 5-bead rod would order at a much higher T* than a tethered nanosphere system where 
the nanoparticle is only a single bead, if the beads interactions were similar. In order to 
compare between the TNS, TNR, and BCP systems we determine the relationship between 
Flory-Huggins interaction parameter, x, and T* for each of the systems. The Flory-Huggins 
interaction parameter allows us to make comparisons between systems regardless of the 
specific model details used. For example, x has been used to compare the calculated phase 
boundaries of BCPs for Brownian dynamics systems that utilize attractive LJ interactions, 
dissipative particle dynamics systems that utilize only repulsive interactions, and mean field 
theory calculations [22]. To determine x, we follow a similar procedure to that outlined in 
References [24J and [22] . For a two-component mixture the free energy can be expressed as: 

~ ln fA + TT ln \X ~ Ia) + ^pr ( 5 ) 



k h T v A N A v B N B {v A v B ) 1/2 

where f A and (1 — f A ) are the fraction of constituents A and B, respectively, v A and v B 
are the volume of the beads of constituents A and B, respectively, and N A and N B are 
the number of beads of constituents A and B, respectively [25]. This expression assumes 
incompressibility where f B = (1 — f A )] note f A should not be confused with the bulk volume 
fraction (j). If we consider the system to be in equilibrium, the free energy will be at a 
minimum and thus = 0- By taking the derivative of the free energy in Equation 5 with 
respect to f A and setting the resulting expression equal to zero, we arrive at the following 
equation that relates x to f A : 

= (vAVB) 1/2 (-ln(f A )v B N B - v_bNb + ln(l - f A )v A N A + v A N A ) 
X v A N A v B N B (2f A -l) [ } 
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A. General procedure for determining \ 

To determine how x scales with T* for each system, we calculate the relative solubility 
of species A mixing into species B as a function of T* and solve for x from Equation [6| 
For example, in the case of the tethered nanosphere system, we specify that the nanosphere 
is species A and the 8 bead polymer chain is species B, and utilize the following general 
procedure. Note that the immiscible components (here, nanospheres and tethers) are not 
bonded to each other for the determination of X- 

Utilizing a long rectangular box with aspect ratio 1:1:4 (x:y:z), we place all the 
nanospheres on one half of the box and all polymers on the other half, creating an interface 
between the two immiscible species. For a specific we run the system for approximately 
10 million timesteps monitoring the potential energy to ensure we reach equilibrium. After 
these 10 million equilibration timesteps we generate an average concentration profile along 
the long dimension of the box by collecting data over the next 10 million time steps. To solve 
for X) we calculate the average fraction of nanoparticles that have mixed into the polymer 
region (i.e. /a) from the concentration profile, then plug this value into Equation |6j solving 
for x- We then repeat for various values of T* creating a relationship between T* and 
X, as shown in Figure [5J Typically, for incompressible mixtures, as T* is decreased x will 
increase and this relationship is often described as x — §^ + D [25] . where C and D are 
system dependent fitting parameters; consequently, we fit our data using linear regression 
to determine the relationship. Since our systems are compressible the x mapping depends 
on the bulk volume fraction, (f). Specifically, the slope of the fitting will increase as (j) is 
increased, i.e. there will be less desire for the two systems to mix. 

B. x mappings for TNS, TNR, and BCP 

Utilizing the same procedure for each of the building blocks, we find the following re- 
lationships for x vs - T* (also shown in Figure [2]). These mappings were performed at 
4>tns = 0.3 and (/)tnr = 0.21, corresponding to the bulk volume fractions where the DG 
was simulated for each of the systems, respectively; the mapping for the BCPs was performed 
at 4>bcp = 0.21 for appropriate comparison with the TNR system. 
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Xtns = (0.93 ± 0.06) /T* - (0.12 ± 0.05) (7) 



XTNR = (6.89 ± 0.53) /T* - (1.28 ± 0.24) (8) 



Xbcp = (6.89 ± 0.15) /T* - (1.70 ± 0.09) (9) 
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FIG. 2: x mappings for TNS (circles), TNR (squares), and BCP (diamonds). 

These relationships follow the expected trends. For an equivalent \ value, T* is higher for 
the TNR system than the BCP system; this corresponds to the behavior seen in reference 
[T2] where the TNR system microphase separated at a higher T* than the equivalent BCP 
system. For an equivalent \ value, T* is substantially higher for the TNR system as 
compared to the TNS system; previous simulations of TNRs found the system microphase 
separates at approximately T*=l [12] whereas the TNS system was found to order at about 
T^=0.3 [H]. We should be aware that these results are for moderate bulk volume fractions 
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(0 = 0.3 for TNS and = 0.21 for TNR and BCP); the results in References [24J and 
[22] show good agreement with theory using this procedure, however the comparisons were 
made at higher values of (j) corresponding to melt conditions (0 = 0.45). The difference in 
(ft will manifest itself in the compressibility assumption, specifically fs = (1 — /a) may not 
be completely valid at the lower (j) values we utilize. However, since our simulations are for 
building blocks in selective solvent (i.e. tether beads are treated with the repulsive WCA 
potential) , the tethers can fill the available space without forming low density holes and thus 
the x mapping should allow for reasonable comparisons among these systems at the lower 
values of (j) . 



III. RESULTS 



In previous publications we reported the presence of the double gyroid (DG) structure in 
both the TNS [Hj and TNR [15] systems. The double gyroid is a bicontinous structure where 
the minority component, in our case the nanoparticles, forms two distinct, interpenetrating 
networks that never connect. The minority component organizes into a series of cylindrical 
tubes (arms) where three tubes connect at each node. Figure [3] shows a simulation snapshot 
of the DG formed by the TNS system (a) and the TNR system (b). For clarity we have 
removed the tethers and show only the minority component (nanospheres or nanorods) 
where the two distinct networks are colored red and white; the particles in the red domain 
are chemically identically to those in the white domain. These structures were identified 
both visually and by calculating the structure factor. The structure factor for the DG shows 
two strong peaks with a characteristic ratio of \/3:Vi as expected [26j. Averaging over 
several runs, the DG was found to form for (j) = 0.3 at 1/T* > 3.28 for the TNS system and 
for (j) = 0.21 at 1/T* > 0.825 for the TNR system. Recent theoretical predictions of the 
order-disorder transition for TNS agree with our simulations, reporting a value of 1/T* ~ 



3 for cj) = 0.3 [27]. Using the x mappings we calculated in section II B, we find reasonable 
agreement between the average order-disorder transitions for both systems as expected. 
Specifically, we find that XorTr-disorder = 3.17 ±0.25 and X™-^Wer = 4.40 ±0.70, where 
the reported error for x is calculated from the maximum error in the linear regression of x 
vs. T* 
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FIG. 3: (a) DG phase formed by the TNS system, (b) DG phase formed by the TNR system. In 
both cases tethers have been removed for clarity. The images correspond to the minimal unit cell 
duplicated once in each direction for clarity. 

A. Packing frustration analysis 

The formation of the DG in the tethered nanoparticle systems is surprising since it is 
known to exist only in a small region of the phase diagram for BCPs [281 [291 EDI EI], rod-coil 
BCPs [16j [32], and surfactants [2], and was not seen in some of our previous simulations of 
the TNS [13] and TNR [12] systems. The limited range of stability of the DG phase in BCP 
systems has been attributed to packing frustration at the nodes of the gyroid [33], EI]- It 
has been shown that the standard deviation in mean curvature, gh , correlates to packing 
frustration and dictates the overall stability of a structure [33] . For example, Matsen and 
Bates [33] calculated gh for various structures finding gh = 0.121 for the gyroid as compared 
to &h = 0.003 for cylinders and attributed this difference to the inability of the gyroid 
structure to simultaneously minimize surface area and minimize packing frustration [33] . 
The magnitude of gh is important in determining which phase will form as the system 
should selectively prefer to form a structure with less packing frustration [33] : this is why 
BCPs typically form the DG over, for instance, the bicontinuous double diamond structure 
where gh = 0.311 [33]. In BCPs, packing frustration has been shown to manifest itself 
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as a high void fraction (low packing density) within the nodes of the gyroid [35] where 
the polymer needs to stretch to fill the volume dictated by the interface [33J. Martinez- 
Veracoechea and Escobedo have shown that packing frustration in the nodes of BCPs can 
be reduced by adding monomer and homopolymer to the system, increasing the stability of 
phases like the DG [35j and stabilizing other structures such as double diamond [36] . 

To assess the stability of the DG structure for the TNS and TNR systems, we examine 
packing frustration by calculating the relative trends in void fraction at the nodes and 
comparing this with the void fraction of the arms, since measuring a H is problematic in a 
system of discrete particles. To look at this trend we approximate the center of a node and 
calculate the void fraction within a spherical volume drawn from the center, repeating for 
various sphere radii, i cut . If we consider a very large value of r cut) we capture the bulk void 
fraction of the system and as i cut is decreased, we get an increasingly localized picture of 
what occurs at the nodes. We perform the same procedure on the arms, again, capturing 
a relative trend in void fraction. For both the TNS and TNR systems, the void fraction 
of the arms and nodes are nearly identical, as shown in Figure |4j suggesting there is a 
uniform density throughout the DG in both structures, irrespective of whether we are at 
the node. Hence, we do not have a characteristic high void fraction within the nodes, as has 
been shown for BCPs [35], thus the nanospheres and nanorods have reduced the packing 
frustration as compared to a flexible BCP. Additionally, we can compare the void fraction 
between the nodes and the bulk system. For both the TNS and TNR systems, we find that 
for large values of i cut the void fraction is approximately that of the bulk system and as 
i cut is decreased, the void fraction is also decreased, as shown in Figure [4} Specifically, for 
the TNS system, the void fraction for large values of r cut is ^0.7 and for small values ^0.5. 
For the TNR system, the void fraction for large values of r cut is ~0.79 and for small values 
^0.55. Thus particles within the nodes pack more densely than the bulk system as a whole. 
Martinez- Veracoechea and Escobedo used a similar analysis to compare a monodisperse BCP 
system to a blend of two different length BCPs [35]. In the blend, the authors found that 
the longer of the two polymers occupied the nodes of the DG, resulting in a larger range 
of stability compared to the monodisperse system [35j. For the monodisperse system, the 
authors found that the void fraction is higher than the bulk for small values of i cut and 
for the blend the authors found that, like our systems, the void fraction is lower than the 
bulk for small values of i cut [35] . This similarity in trends suggests that the DG structures 
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formed by the TNS and TNR systems may be more stable than an equivalent flexible BCP 
system. We conclude that by not forming low density regions at the nodes of the DG both 
the TNS and TNR systems have reduced their packing frustration, and we hypothesize that 
this reduction is a direct result of the geometry of the minority component and its ability to 
pack locally into compact, low energy structures. We test this hypothesis in sections [III B 
and UTTCl 



|0.7 
Jo.6 

?0.5 

> 



§0.8 
|0.7 

5 0.6 

o 
> 



(a) 



(b) 



1,1,1,1,1,1,1,1 


i 1 i 1 i 1 i 1 i 1 i 


■ ■ ■ Bulk 

oo Node 

i i i 


i 



5 6 7 8 9 

•cut (°) 



5 6 7 8 9 

r cu. (°) 



1 1 1 1 1 1 1 1 1 1 1 













••• Bulk 
■-■ Arm 






oo Node 




1,1,1,1,1, 


i , i 





FIG. 4: Relative void fraction within the nodes of the (a) TNS DG for 1/T* =3.3 (x = 3.19 ± 
0.25) and (b) TNR DG for 1/T* =0.9 (x = 4.92 ± 0.72). 



We can further assess the importance of the particle geometry on stability by specifically 
comparing to a BCP system. If we replace the rigid constraints in the TNR system with 
FENE springs, we arrive at a simple model for a flexible coil-coil BCP, as was previously 



described in section |I C[ If the rigidity of the nanorod does not influence the stability and 
packing frustration of the DG phase, we would expect to find the DG in the BCP system 
at x ^ 4.40 ± 0.70 and (j) = 0.21. Starting from the ordered DG configuration, we run 
the system as a coil-coil BCP for various values of x- We find that the DG structure 



13 



is not stable for values of x <^ 8; instead, the DG falls apart and forms a disordered 
aggregate. An example of the disordered structure is shown in Figure [5j For values of 
X >~ 8, the DG structure persists, however, it is most likely a kinetically arrested structure 
and not at equilibrium due to the large value of x- Starting from an athermal, disordered 
configuration, we incrementally cool the system, finding only disordered structures, even for 
values of x >^ 8. This supports the contention that the DG is not stable for the BCP 
system under these conditions. Note that simulations were also performed with various unit 
cell sizes to avoid any box size issues that are associated with 3d periodic microstructures, 
again, yielding only disordered structures. Thus, we find that for equivalent statepoints, we 
were unable to realize the DG phase when the rigid constraint was removed, which suggests 
that the stabilization of the DG is strongly influenced and controlled by the geometry of the 
aggregating species. 




FIG. 5: Disordered structure formed by BCP system for x — 8-6 ± 0.26. The tethers have been 
removed for clarity, and only the minority, aggregating species is present. 

B. Local structure of the TNS gyroid 

The ability of the TNS system to reduce the packing frustration in the DG can be un- 
derstood by looking at the local packing of the particles at the node. Figure [6] shows a 
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FIG. 6: Node of the DG formed by the TNS system, with a perfect icosahedron showed in inset. 
Tethers have been removed for clarity. 

simulation snapshot of the TNS DG node. We see distinct local ordering in the node, specif- 
ically ring-like structures that resemble icosahedral clusters. An icosahedron is constructed 
of a central particle surrounded by 12 nearest neighbors and is a minimal potential energy 
structure for 13 Lennard- Jones particles. This ordering is not limited to the nodes and 
occurs throughout the entire DG structure. 

In reference [2] we introduced the R y i m method based on spherical harmonics [37J and 
used this to identify the local packing of these particles. The general R y i m method relies on 
a rotationally invariant spherical harmonic "fingerprint" of the structure. In this method 
we first calculate q™ by summing over all nearest neighbor directions between particles: 



where fij is the vector drawn from particle i to its nearest neighbor j, is the total number 
of neighbors, £ is the specific harmonic, and Y^ m is the spherical harmonic expansion [37] . 
From this we can construct two rotationally invariant measures of cluster shape, Q^, the 
vector magnitude of qj 1 , and W£ : the rotationally invariant combination of the average values 
of q™ [37] . As defined by Steinhardt, et al : 
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where [ is the Wigner 3j symbol [37]. 

mi m 2 m 3 / 

To determine the local structure of a particle we calculate the Qi and W£ values for / = 
4, 6,... 12, and then calculate the residual value with respect to a library of known structures 
with matching coordination {H]: 



R i = \IYaIa (Ql ~ Qrefif + E?=4 ( W l ~ W reff • 



A particle is considered to be in the local configuration i that minimizes the residual Ri or 
considered to be disordered if the residual value exceeds a certain cutoff, chosen as 0.316 for 
this specific system and conditions. The R y i m method is well suited to accurately identify 
particle configurations as it provides a rotationally invariant description of a local structure 
and does not rely on multiple cutoff values to determine the configuration, minimizing 
potential error. Additional description of the method can be found in references [T4l [38] . 




cn = 12 cn = 11 cn = 10 cn = 9 cn = 8 

FIG. 7: Icosahedral clusters ranging from full coordination number (cn) of 12 to partial coordination 
of 8. 

Our reference library includes standard particle arrangements such as the Kasper poly- 
hedra (Z10, Zll, Z12, Z13, Z14, Z15) [39J, face centered cubic (FCC), hexagonally close 
packed (HCP), body centered cubic (BCC), and simple cubic (SC). Additionally, our li- 
brary includes partial icosahedral clusters, where 1-4 particles are removed from the ideal 
Z12 Icosahedron (see Figure [7]), and partial clusters of FCC and HCP where 1-5 particles 
are removed from the ideal clusters. 
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Using the R y im method, we determined in reference [H] the packing of the nanoparticles 
in the DG to be predominantly icosahedral clusters with partial coordination, i.e. we found 
clusters that retain the same bond angles as a perfect icosahedron but with 1-4 particles 
removed (see Figure [7| [14] ; these clusters are identical to the minimal potential energy clus- 
ters found by Doye and Wales [40] . Partial clusters are formed as a result of the steric effects 
of the tethers and the microphase separation that occurs. Table [I] shows the percentage of 
nanospheres that are "central" particles within partial icosahedral clusters, the percentage 
of nanospheres that are "central" particles within partial crystalline cluster (HCP and FCC 
in this case) along with the average coordination number of nanoparticles for an example 
cooling sequence. The double line in Table |T] signifies the transition from a disordered state 
to the ordered DG microphase. As we increase l/T* (i.e. cool the system), we notice 
a substantial increase in the number of icosahedral-like clusters, increasing from approxi- 
mately 17.5% to 30% but very little change in crystalline arrangements, with an average 
value approximately 16% within the ordered regime (l/T* > 3.25). We also see a minor 
increase in average coordination number, increasing from approximately 7 to 8; this does 
not change as rapidly as the percentage of icosahedral clusters since coordination number 
does not differentiate between ordered and disordered local configurations. 



TABLE I: Local Structure of Nanospheres 



l/T* 


X ± error 


% Icos. ± stdev 


% Crystal ± stdev 


Coord. Num. ± stdev 


3.00 


2.91 ± 0.23 


17.3 ± 1.97 


17.1 ± 1.76 


6.9 ± 2.3 


3.15 


3.05 ± 0.24 


18.5 ± 2.91 


18.1 ± 2.53 


7.2 ± 2.1 


3.25 


3.14 ± 0.25 


21.7 ± 2.67 


18.7 ± 2.66 


7.4 ± 2.0 


3.30 


3.19 ± 0.25 


23.9 ± 1.81 


16.9 ± 1.91 


7.4 ± 2.1 


3.50 


3.38 ± 0.26 


26.4 ± 2.80 


16.2 ± 3.13 


7.5 ± 2.0 


3.60 


3.47 ± 0.27 


28.7 ± 2.01 


15.7 ± 2.49 


7.6 ± 2.0 


3.70 


3.56 ± 0.27 


28.6 ± 3.12 


15.8 ± 2.34 


7.6 ± 2.1 


3.80 


3.66 ± 0.28 


30.1 ± 2.61 


15.5 ± 1.62 


7.8 ± 2.0 



In Figure [8] we also present these results grouped by the coordination of the partial 
icosahedral cluster, i.e. we plot the data for each cluster arrangement separately. As the co- 
ordination number suggests, at all T* we are most likely to find clusters with a coordination 
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number of 8 and likewise, we find only a very small number of particles with a coordination 
number of 11. The summation of all of the icosahedral clusters, as reported in Table []] is 
shown as solid diamonds in Figure [8] demonstrating a clear linear increase over the range 
considered. 




FIG. 8: Icosahedral clusters with partial coordination. Data is grouped by the coordination number 
of the cluster, ranging from 8 to 11. The solid diamonds correspond to the sum of all partial 
icosahedral coordinations, also reported in Table [TJ All data was fit using a linear regression. The 
error bars correspond to the standard deviation. 

Similarly, we can group the results for crystalline partial clusters by their coordination 
number. This eliminates some complexity, since for HCP and FCC there are 12 unique 
partial clusters for coordination numbers ranging from 7-11; this is due to the fact that the 
partial cluster we form depends on which particle(s) we choose to remove. For example, 
we have two unique spherical harmonic fingerprints for an HCP cluster with coordination 
number of 10 depending on which particles we remove. Additionally, because of the similar- 
ities between HCP and FCC, certain partial clusters are identical and we cannot determine 
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whether the cluster is FCC or HCP, only that it possesses characteristics of both. As such, 
we group together the partial clusters of both FCC and HCP by coordination number of the 
partial cluster, as shown in Figure [9j We find that for all configurations, the percentage of 
crystalline clusters does not change much as we increase 1/T* (i.e. cool the system); there 
is a minor decrease in the percentage of clusters on the order of the standard deviations. 
We also see that we are most likely to have crystalline partial clusters with coordination 
numbers of 9 and 7; as we saw for icosahedral clusters, we find only a small amount of 
clusters with large coordination numbers of 11. 




FIG. 9: Percentage of nanospheres that are central particles in FCC and HCP clusters with partial 
coordination. Data is grouped by the coordination number of the cluster, ranging from 7 to 11. 
The solid diamonds correspond to the sum of all partial crystal coordinations, also reported in [IJ 
All data was fit using linear regression. The error bars correspond to the standard deviation. 

The fact that icosahedral arrangements are favored over crystalline ordering is somewhat 
surprising, as we find for a large bulk system of nanopsheres without tethers that FCC/HCP 
crystalline ordering is dominant; icosahedral arrangements are favored only for simulations of 
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small numbers of particles. In reference [14J we found that icosahedral ordering of nanopar- 
ticles was a result of the confinement that occurs as the system phase separates. We showed 
that by confining particles into cylinders with hard walls, there was a transition from pre- 
dominantly crystalline ordering to predominantly icosahedral ordering when the diameter 
of the cylinder (scaled by particle size) was less than 5, corresponding to the approximate 
diameter of the domains in the DG system [14J. 

C. Local structure of the TNR gyroid 



FIG. 10: Node of the DG formed by the TNR system with a hexagonal bundle highlighted. Tethers 
have been removed for clarity. 



Similar to the TNS system, the DG formed by the TNR system also has distinct local 
ordering of the nanoparticles in addition to the bulk microphase separation. Figure [TO] shows 
a simulation snapshot of the TNR node. We notice that the nanorods attempt to adopt 
bundled structures, where a nanorod is surrounded by six nearest neighbors in a hexagonal 
fashion; a single bundle is highlighted in Figure [To} A six neighbor bundle is the densest, 
minimal potential energy structure for seven rods, representing full coordination, analogous 
to a coordination of 12 for an icosahedron. The tendency to form these bundles can be 
observed by examining the histogram of coordination number of the center-of-mass of each 
rod, shown in Figure 11 for 1/T* = 0.9. There is a clear bias towards high coordination 



numbers and we find no coordinations greater than six. We do find partially coordinated 
clusters as a result of rods being situated on the boundary with the tether region, as was 
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also the case of the TNS system. Example clusters from our simulations are shown in 



Figure [12] for coordination numbers, cn = 3 - 6, where the preference to hexagonal packing 
is highlighted. 

To reduce the grafting density of the tethers (i.e. the local density of tethers in a small 
region), thus maximizing entropy for the tether, we expect the grafting points of the tethers 
(i.e. the points in 3d space where the tether is attached to the nanorod) to be equally 
distributed along the interface between the nanorods and tethers. In Figure [TT] we plot the 
histogram of the coordination number of the grafting points noticing a strong tendency for 
coordinations of one and two and do not find any coordinations above four. A bundle of 
seven rods would have the most unfavorable configuration in terms of entropy if the grafting 
point coordination number were six, corresponding to all tethers being oriented on the same 
side of the bundle. The histogram shows that each bundle has two or three tethers oriented 
in the same direction per bundle and thus tether attachment is well distributed. 

We observe there is a deviation from the ideal hexagonally packed bundle structure as a 
result of the tether attachment; the bundles tilt and splay with respect to the central nanorod 
additionally reducing the grafting density of the tethers [12] ; this behavior is evident in the 



clusters in Figure 12, This manifests itself as a twisted structure in the tubes (arms) of 
the DG, as has been seen in previous simulations of tethered nanorods that form twisted 
cylinders [12]. To quantify this behavior, we calculate the angle between the director of a 
nanorod with that of its nearest neighbors by calculating the dot product. On average we 
find that the angle between a nanorod and its neighbors is 10.7 degrees with a standard 
deviation of 6.1 degrees for a system with 1/T* = 0.9. Alternatively, we can quantify the 
splay of the rod bundles by calculating the nematic order parameter, S, for each bundle. To 
calculate the nematic order parameter, we first calculate the 3x3 nematic order tensor, 

N 



(?«-¥) 

i=l v 7 



where a, (3 = x,y,z, u l a is the a component of the rod director and 5 a p is the Kronecker delta 
[4T] . We take S to be the largest eigenvalue of the matrix Q [41]. This construction results 
in perfectly crystalline systems having a value of 5=1, nematic ordered liquid crystalline 
systems having S—0. 3-0. 9, and isotropic systems having S<0.3. We find that for a system 
with 1/T* — 0.9, the bundles of nanorods have an average value of S— 0.96, demonstrating 
a small deviation from the ideal crystalline behavior. Note that this is the average local 
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FIG. 11: Histogram of the coordination number of centers-of-mass of the rods (grey) and a his- 
togram of the coordination number of grafting points (black) for 1/T* = 0.9. 




cn = 6 cn = 5 cn = 4 cn = 3 



FIG. 12: Example clusters formed by rods at 1/T* = 0.9 for coordination numbers 6, 5, 4, and 3. 

nematic order parameter of the individual bundles, not the order parameter of the bulk 
system. An alternative form for the nematic order parameter is given by S = (l/2)(3cos 2 9 — 
1), where 9 is the angle between a given nanrod and the director [25J. Substituting the 
average splay angle of 10.7 degrees into this form yields a value of 5=0.95 showing good 
agreement between these two methods of quantifying the splay. 



In Table II we calculate both the average splay angle and nematic order parameter of 
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the bundles for an example cooling sequence, where the system is equilibrated at each 
1/T* prior to subsequent cooling; the transition from disordered to DG is denoted by the 
double line. We notice that for the ordered microphases (1/T* > 0.8) the splay angle is 
roughly 10 degrees, with a standard deviation of 7.2 degres for the lowest value of 1/T*, 
and 3.8 degrees for the highest value of 1/T*. Within the ordered region, the nematic order 
parameter deviates little from the average value of 0.965. Within the disordered region 
(1/T* < 0.8), both the splay angle and standard deviation are large, suggesting that there 
is a large variation in alignment of a set of neighboring rods; this is supported by values of S 
that are clearly less than the average value for the ordered DG configurations, indicating that 
neighboring rods are not strongly aligned. The combination of splay angle and S allows us to 
conclude that, similar to the TNS system, there is a strong connection between microphase 
separation and local packing of the nanoparticles; we find strong local ordering where we 
find the DG, and very little local ordering elsewhere. 



TABLE II: Local coordination in TNR system 



1/T* 


X ± error 


Splay Angle (degrees) ± stdev 


S 


0.50 


2.17 ± 0.51 


44.7 ± 25.7 


0.73 


0.70 


3.54 ± 0.61 


34.8 ± 25.7 


0.79 


0.80 


4.23 ± 0.66 


19.6 ± 17.3 


0.88 


0.85 


4.58 ± 0.70 


11.4 ± 7.2 


0.96 


0.90 


4.92 ± 0.72 


10.7 ± 6.1 


0.96 


0.95 


5.27 ± 0.75 


10.0 ± 4.7 


0.97 


1.00 


5.61 ± 0.77 


10.0 ± 4.7 


0.97 


1.10 


6.30 ± 0.83 


9.7 ± 4.9 


0.97 


1.20 


6.99 ± 0.88 


9.4 ± 3.8 


0.97 



We also calculate the rotation of the rod as a function of time and 1/T* to assess the 
ability of the rods to reorient themselves. We calculate the average rotation of the rods by 
calculating the dot product of a rod's director at the current time with its initial starting 



director, as plotted in Figure 13, We find that at small values of 1/T* there is a large 
amount of rotational freedom in the rods; we note that the rods for 1/T* =0.5 have, on 
average, rotated approximately 85 degrees at time = 100 (10000 timesteps). The large 
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amount of rotational mobility explains the large splay angle and large standard deviation 
for high temperature disordered phases, as reported in Table |TI| This rotational freedom is 
decreased as we increase 1/T* (i.e. cool the system); it is clear that for values of 1/T* > 
0.8, where the system has microphase separated into the DG, there is a drastic reduction 
in this rotational freedom and the average rotation of rods in the system is less than 15 
degrees over the entire time range sampled. The lack of rotational mobility suggests that 
the direction of the rod is strongly correlated through time and, again, this supports the 
results in Table [TT] where, for the DG, the distribution of splay angle is relatively narrow as 
compared to the disordered structure. 



i i i i i i i i i I i i i i i i i i i I i i i i i i i i i I i i i i i i i i i I i i i i i i i i i 
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FIG. 13: Average rotation as a function of time for various values of 1/T*. Solid black circles 
represent disordered microstructure, open circles correspond to the DG morphology. 



IV. CONCLUSIONS 



We have performed Brownian dynamics simulations of tethered nanospheres and tethered 
nanorods that predict the double gyroid structure. In both the tethered nanosphere and 
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tethered nanorod systems we see a common trend towards the formation of dense, minimal 
potential energy structures. Tethered nanospheres predominantly form icosahedral clusters 
with partial coordination. Tethered nanorods form hexagonally packed bundles of rods, 
with both full and partial coordination, where the rods tilt and splay with respect to their 
neighbors. The icosahedral and hexagonal bundle structures adopted by the nanoparticles 
allow for a uniform density throughout the double gyroid structure; we do not see low 
density regions at the nodes of the double gyroid as would be expected for flexible chains 
[35] . The result of these local configurations is a reduction in packing frustration, which has 
been linked to the overall stability of the phase [33j CSJ [36] . We have seen that the rigidity 
of the nanorods is crucial for the stabilization of the double gyroid relative to BCPs; for 
equivalent conditions, we did not find the double gyroid structure for an analogous flexible 
diblock copolymer system. This suggests that by properly choosing the geometry of the 
nanoparticle, we may be able to not only stabilize the double gyroid, but potentially other 
bicontinuous structures by harnessing the unique combination of microphase separation and 
local packing. 
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